%% Solvers for ERME Inner Iterations
% The power method for ERME's consists of a two-level iteration.  The inner
% iteration is an eigenproblem of the form
%
% $$MR(k_m)*J = \lambda*J$$
%
% This can be solved via power iteration (PI):
%
% $$    J_{n}    = MR(k_m)*J_{n-1}$$
%
% $$    \lambda  = || J_{n} ||_2    $$
%
% $$    J_{n}    = J_{n} / \lambda  $$
% 
% While PI works, it's not particular efficient, especially for problems with a
% high dominance ratio.  Here, even for a relatively small model, the dominance
% ratio of MR is roughly 0.994.  Hence, we seek alternative schemes to
% accelerate these inner iterations. To this end, we compare:
%
%  * power method (PI)
%  * successive over-relaxation (SOR) acceleration of PI
%  * SOR with omega via Aitken $\delta^2$ acceleration of PI
%  * Chebyshev acceleration of PI
%  * Implicitly restarted Arnoldi
%  * MATLAB's built-in sparse eigensolver "eigs" for just the dominant eigenpair
%
% The model problem is the two-group IAEA benchmark.  Two evaluations of MR(k) 
% were exported from SERMENT.  The first is evaluated after one outer iteration,
% while the second is the final evaluation for converged k.  The model used
% 17x17 assemblies with 4th order spatial expansions.  Thus, there are
%
% $$ 2 [g] \times (4+1) [orders] \times 17^2 [mesh] \times 4 [faces] = 11560 \: \mathrm{unknowns} \, ,$$  
%
% likely a sufficient size for evaluating algorithms.

clear%, clc
% load MRb;     % MR is first, MRb is last
% A       = MR; % Rename to A for simplicity,
% clear MR;     % and clear MR from memory.
load MRsmall2;
A = MR;
N       = length(A);            % Size of the problem.
[vt,et] = eigs(A,3);            % Get the first three eigenpairs.
                                % L_2 / L_1 = rho, L_3/L_1 = "b"
et      = diag(et);             % put eigval's in a vector
[e,idx] = sort(et,1,'descend'); % sort them, high to low, and save index 
v       = vt; 
for i = 1:length(e)
    v(:,i) = vt(:, idx(i));     % sort the vectors
end

% Reference values.  
L_ref   = e(1);
X_ref   = v(:,1);
rho_ref = e(2) / e(1); % reference dominant ratio and
b_ref   = e(3) / e(1); % ratio of largest negative value to dominant value
                       % which is needed for Chebyshev

% Convergence criteria.                       
maxit = 1e5; 
tol   = 1e-3;
Xo = ones(length(A), 1);  % Set just the zeroth order
Xo(2:2:end) = 0;          % moments to 1
Xo = Xo / norm(Xo);       % and normalize to unity.

%-------------------------------------------------------------------------------

% Power Method
tic, disp(' *** PI *** ')
[X1, L1, X_err1, L_err1, iter1 flag] = ...
    powermethod(A, Xo, X_ref, L_ref, maxit, tol);
fprintf(' '),toc,disp(' '),disp(' ')

% % PI+SOR
% tic, disp(' *** PI+SOR *** ')
% [X2, L2, X_err2, L_err2, iter2, flag] = ...
%     sorpower(A, Xo, X_ref, L_ref, maxit, tol, 1, 1.0003);
% disp(sprintf(' Percent decrease in number of iterations: %6.5f ', ...
%     100*(iter1-iter2)/iter1)),fprintf(' '),toc,disp(' '),disp(' ')
% 
% % Aitken-Accelerated Power Method
% tic, disp(' *** Aitken *** ')
% [X3, L3, X_err3, L_err3, iter3, flag] = ...
%     aitkenpower(A, Xo, X_ref, L_ref, maxit, tol);
% disp(sprintf(' Percent decrease in number of iterations: %6.5f ', ...
%     100*(iter1-iter3)/iter1)),fprintf(' '),toc,disp(' '),disp(' ')

% Chebyshev-Accelerated Power Method
tic, disp(' *** Chebyshev *** ')
[X4, L4, X_err4, L_err4, iter4, flag] = ...
    cheby(A, Xo, X_ref, L_ref, 1e5, tol);
disp(sprintf(' Percent decrease in number of iterations: %6.5f ', ...
    100*(iter1-iter4)/iter1)),fprintf(' '),toc,disp(' '),disp(' ')

%Implicitly Restarted Arnoldi Method
tic, disp(' *** IRAM *** ')
[X5, L5, X_err5, L_err5, iter5, flag] = ...
    arnoldi(A, Xo, X_ref, L_ref, maxit, tol);
toc,fprintf(' '),disp(' '),disp(' ')


%Inverse Iteration 
% tic, disp(' *** Inverse *** ')
% [X6, L6, X_err6, L_err6, iter6, flag] = ...
%     invit(A, Xo, X_ref, L_ref, maxit, tol);
% disp(sprintf(' Percent decrease in number of iterations: %6.5f ', ...
%     100*(iter1-iter6)/iter1)),fprintf(' '),toc,disp(' '),disp(' ')


% % Matlab "eigs" (probably IRAM variety)
% tic, disp(' *** eigs *** ')
% [X6, L6] = eigs(A, 1);
% toc
% % 


